Setup

Set up workspace, load data, and load required packages.

rm(list=ls(all=TRUE))

results <- read.csv("Data/Urchin_ClimateChange_Data.csv", header=T, na.strings="NA")

library("lme4") #linear mixed modeling
library("lmtest") #linear mixed modeling
library("lmerTest") #calculate p-values in models
library("car") #levenes tests
library("emmeans") #post-hoc tests
library("lsmeans") #post-hoc tests
library("plotrix") #calculate std.error
library("plyr")  #splitting, applying, and combining data
library("dplyr") 
library("cowplot") #grid plotting and arrange plots
library("effects") #plot effects of modeling
library("MuMIn")
library("car") #levenes tests, Anova, qqPlots, 
library("tidyverse")
library("stats")
library("onewaytests") #allows for Welch test with unequal variances (spine lengths)
library("stats")
library("effsize")
library("FSA")
library("rcompanion")
library("broom") #turns model output into data frames
library("kableExtra") #makes data tables easier to read
library("tidyverse")
library("pwr") #run power tests

Approach:

1) Build and run model  
2) Check for normality of residuals  
3) Check for homogeniety of variance of residuals  
4) Look at model summary  
5) Run anova to check for significance of main effects  

Analysis:

  1. Environmental Conditions and Manipulations:
    1. Calculate means
    2. Check environmental manipulations with ANOVAs for treatments
    3. Figure 1: Temp over time
    4. Figure 2: pH over time
  2. Urchin Diameters:
    1. Calculate diameter means for:
      1. Day -24
      2. Day -13 iii.Day 1
      3. Day 126
      4. Table 3: Diameter Means
    2. Check for differences in urchin diameters from different tanks on Day 1
  3. Response Analyses:
    1. Growth (%)
      1. Growth Analysis
      2. Figure 3: Growth Interaction Plot
    2. Calcification Ratio
      1. Calcification Analysis at the Tips
      2. Calcification Analysis at the Bases
      3. Figure 4: Calcification Interaction Plot
    3. Spine Loss (count)
      1. Spine Loss Analysis
      2. Figure 5: Spine Loss Interaction Plot

1. Environmental Conditions and Manipulations:

a. Calculate environmental means

Table 1: Environmental Means
EnviroMean <- 
  results %>% 
    select("Day", "Temperature", "pH", "Temp","pHspec", "pCO2out", "AlkTotal") %>%
    filter(Day >= "1", Day <= "123") %>%
    drop_na(Day) %>% 
    group_by(Temperature, pH) %>% 
    summarise(mean_T = mean(Temp, na.rm = T),
              se_T = sd(Temp, na.rm = T)/sqrt(6),
              min_T = min(Temp, na.rm = T), 
              max_T = max(Temp, na.rm = T),
              mean_pH = mean(pHspec, na.rm = T), 
              se_pH = sd(pHspec, na.rm = T)/sqrt(6),
              mean_pCO2 = mean(pCO2out, na.rm = T), 
              se_pCO2 = sd(pCO2out, na.rm = T)/sqrt(6),
              mean_AlkTotal = mean(AlkTotal, na.rm = T),
              se_AlkTotal = sd(AlkTotal, na.rm = T)/sqrt(6))

EnviroMean %>% 
  kbl(caption = "Table 1: Experimental treatments and environmental conditions (mean +/- standard error) of the 126-day experiment.") %>% 
  kable_classic_2() %>% 
  kable_styling(c("striped", "hover"), full_width = F, html_font= "Times New Roman")
Table 1: Experimental treatments and environmental conditions (mean +/- standard error) of the 126-day experiment.
Temperature pH mean_T se_T min_T max_T mean_pH se_pH mean_pCO2 se_pCO2 mean_AlkTotal se_AlkTotal
ambient acidified 24.79000 0.8869310 22.07 27.36 7.644120 0.0101512 1082.9250 65.31879 2142.000 112.494347
ambient ambient 24.78429 0.8856820 22.09 27.40 7.961042 0.0146725 480.8325 22.56737 2181.872 5.679540
heated acidified 26.73429 0.8269073 23.75 29.38 7.662486 0.0129196 1106.5750 77.15382 2137.082 123.057409
heated ambient 26.70434 0.8184090 23.82 29.34 7.956690 0.0141496 524.3757 20.91815 2180.797 5.838938

b. Check environmental manipulations with ANOVAs

Environmental <- 
  results %>% 
    select("Day", "Temp", "pH", "Temperature","pHspec","Header", "TankID") %>%
    filter(Day >= "1", Day != "126") #use only days from start of experiment when desired conditions were reached
  

#temp between high and ambient temp
tempmod<-aov(Temp ~ Temperature, data=Environmental)
summary(tempmod) #temperature was different between high and ambient treatment

#pH between high and ambient co2
summary(aov(pHspec ~ pH, data=Environmental)) #pH was different between high and ambient co2 treatments



##HEADERS
#temp between headers in high
aov1<- Environmental %>%
  filter(Temperature=="heated")

modaov1<-aov(Temp ~ as.factor(Header), data=aov1)

summary(modaov1) #not different between headers

#temp between headers in ambient temp
aov2<- Environmental %>%
  filter(Temperature=="ambient")

summary(aov(Temp ~ as.factor(Header), data=aov2)) #not different between headers

#pH between headers in acidified
aov3<- Environmental %>%
  filter(pH=="acidified")

summary(aov(pHspec ~ as.factor(Header), data=aov3)) #not different between headers

#pH between headers in ambient pH
aov4<- Environmental %>%
  filter(pH=="ambient")

summary(aov(pHspec ~ as.factor(Header), data=aov4)) #not different between headers


##TANKS
#temp between tanks in high temp
aov1<- Environmental %>%
  filter(Temperature=="heated")

summary(aov(Temp ~ as.factor(TankID), data=aov1)) #not different between heated tanks

#temp between tanks in ambient temp
aov2<- Environmental %>%
  filter(Temperature=="ambient")

summary(aov(Temp ~ as.factor(TankID), data=aov2)) #not different between ambient tanks

#pH between tanks in acidified
aov3<- Environmental %>%
  filter(pH=="acidified")

summary(aov(pHspec ~ as.factor(TankID), data=aov3)) #not different between tanks in acidified

#pH between headers in ambient pH
aov4<- Environmental %>%
  filter(pH=="ambient")

summary(aov(pHspec ~ as.factor(TankID), data=aov4)) #not different between headers
Table 2: Results of ANOVA for environmental
Treatment Comparison Pr(>F)
Heated v. Ambient <0.001
Acidified v. Ambient <0.001
Ambient temp headers 0.998
Heated headers 0.842
Ambient pH headers 0.129
Acidified headers 0.150
Ambient temp tanks 1.000
Heated tanks 0.999
Ambient pH tanks 0.675
Acidified tanks 0.888

Check assumptions of treatment conditions model

#check assumptions
# 1. Normality of residuals
qqPlot(residuals(tempmod)) #very much not normal

shapiro.test(residuals(tempmod)) # nope
hist(residuals(tempmod)) #nope...

#check assumptions
# 2. Equal variances
#leveneTest(residuals()~$)
#plot(fitted(),resid(,type="pearson"),col="blue")

Check using nonparametric Kruskal-Wallis test

# did i set this up right?
kruskal.test(Environmental$Temp~Environmental$Temperature)# **diff
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Environmental$Temp by Environmental$Temperature
## Kruskal-Wallis chi-squared = 203.29, df = 1, p-value < 2.2e-16
kruskal.test(Environmental$pHspec~Environmental$pH)# **diff
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Environmental$pHspec by Environmental$pH
## Kruskal-Wallis chi-squared = 470.87, df = 1, p-value < 2.2e-16
kruskal.test(Environmental$pHspec~Environmental$Temperature)# Not dif
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Environmental$pHspec by Environmental$Temperature
## Kruskal-Wallis chi-squared = 0.080493, df = 1, p-value = 0.7766
kruskal.test(Environmental$Temp~Environmental$pH)# not dif
## 
##  Kruskal-Wallis rank sum test
## 
## data:  Environmental$Temp by Environmental$pH
## Kruskal-Wallis chi-squared = 0.18832, df = 1, p-value = 0.6643

c. Figure 1: Temperature over time

tempsummary<-results %>% 
  select("Day", "Temperature", "Temp") %>%
  drop_na(Temp) %>%
  group_by(Day, Temperature) %>%
  mutate(meanT = mean(Temp)) %>%
  group_by(Temperature) %>% 
  mutate(sd = sd(Temp)) %>%
  mutate(se = sd/sqrt(12)) 

tempplot<-ggplot(data=tempsummary, aes(Day, meanT, color = Temperature)) +
    geom_point(size = 2.5, show.legend=FALSE) +
    geom_errorbar(aes(ymin = meanT-se, ymax = meanT+se)) +
    geom_line(size=1.2) +
    geom_vline(xintercept=0) +
    scale_x_continuous(name="Time (Day)", breaks = seq(0,130, by = 10)) +
    scale_y_continuous(name="Temperature (°C)", breaks = seq(20,30, by = 1)) +
    scale_shape_discrete(name = NULL,
                         labels = c("Ambient", "High"))+
    scale_color_manual(name = NULL,
                       values = c("gray40", "firebrick3"),
                       labels = c("Ambient", "High")) +
    ggtitle("")+
    theme_classic() +
    theme(plot.margin=unit(c(0.3,0.6,0.3,0.3), "cm"))+
    theme(legend.margin = margin(0), 
          legend.position = c(.98,.8), 
          legend.justification = c("right", "top"),
          legend.background = element_blank(), 
          axis.title = element_text(size = 16, face="bold"),
          axis.text = element_text(size = 14), 
          plot.title = element_text(size=16, face="bold"),
          legend.text = element_text(size=14, face="bold"));tempplot

#ggsave(filename="Figures/20200611/TempFig1.png", plot=tempplot, dpi=300, width=7, height=5, units="in")

d. Figure 2: pH over time

pHsummary<-results %>% 
  select("Day","pH", "pHspec") %>%
  drop_na(pHspec) %>%
  group_by(Day, pH) %>%
  mutate(meanpH = mean(pHspec)) %>%
  group_by(pH, meanpH) %>% 
  mutate(sd = sd(pHspec)) %>%
  mutate(se = sd/sqrt(12)) 
  
  
pHplot<-ggplot(data=pHsummary, aes(Day, meanpH, color = pH)) +
    geom_point(size = 2.5, show.legend=FALSE) +
    geom_errorbar(aes(ymin = meanpH-sd, ymax = meanpH+sd)) +
    geom_line(size=1.2) +
    geom_vline(xintercept=0) +
    scale_x_continuous(name="Time (Day)", breaks = seq(0,130, by = 10)) +
    scale_y_continuous(name="pH (Total Scale)", breaks = seq(7.5,8.1, by = .1)) +
    scale_shape_discrete(name = NULL,
                         labels = c("Low pH", "High pH"))+
    scale_color_manual(name = NULL,
                       values = c("skyblue3", "gray40"),
                       labels = c("Acidified", "Ambient"),
                       guide = guide_legend(reverse = TRUE)) +
    ggtitle("")+
    theme_classic() +
    theme(plot.margin=unit(c(0.3,0.6,0.3,0.3), "cm"))+
    theme(legend.margin = margin(0), 
          legend.position = c(.98,.8), 
          legend.justification = c("right", "top"),
          legend.background = element_blank(), 
          axis.title = element_text(size = 16, face="bold"),
          axis.text = element_text(size = 14), 
          plot.title = element_text(size=16, face="bold"),
          legend.text = element_text(size=14, face="bold"));pHplot

#ggsave(filename="Figures/20200611/pHFig1.png", plot=pHplot, dpi=300, width=7, height=5, units="in")

2. Urchin Diameters:

a. Calculate diameter means

i. Day -24

Assemble data for diameters of initial collection of urchins from hatchery and calculate mean

##Initial Collection = Day -24
InitialDiameter <-
  #seperate out the initial sizes of urchins on day -24 after initial collection.
  results %>% 
    select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
    filter(TankID != "8t" ) %>% 
      #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
    filter(Day == "-24") %>% 
      #sizes on day -24  
    select("TankID", "Diameter", "Temperature", "pH", "Treatment")


InitialMean <- mean(InitialDiameter$Diameter)
InitialSE <- sd(InitialDiameter$Diameter)/sqrt(23)

InitialMean
InitialSE

ii. Day -13

Assemble Data for diameters of acclimated urchins before conditioning and calculate mean

##After acclimating before ramp up = Day-13

AcclimatedDiameter <-
  #seperate out the initial sizes of urchins on day -24, when desired conditions were met.
  results %>% 
    select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
    filter(TankID != "8t" ) %>% 
      #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
    filter(Day == "-13") %>% 
      #sizes on day -24  
    select("TankID", "Diameter", "Temperature", "pH", "Treatment")


AcclimatedMean <- mean(AcclimatedDiameter$Diameter)
AcclimatedSE <- sd(AcclimatedDiameter$Diameter)/sqrt(23) 

AcclimatedMean
AcclimatedSE

iii. Day 1

Assemble Data for diameters on day 1 of experiment when desired conditions were reached

##After ramp up and desired conditions are reached to begin experiment = Day 1

Day1Diameter <-
  #seperate out the initial sizes of urchins on day 1, when desired conditions were met.
  results %>% 
    select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
    filter(TankID != "8t" ) %>% 
      #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
    filter(Day == "1") %>% 
      #sizes on day 1  
    select("TankID", "Diameter", "Temperature", "pH", "Treatment")


Day1Mean <- mean(Day1Diameter$Diameter)
Day1SE <- sd(Day1Diameter$Diameter)/sqrt(23)

Day1Mean
Day1SE

iv. Day 126

Assemble Data for diameters on day 126 of experiment after full growth

##After ramp up and desired conditions are reached to begin experiment = Day 1

Day126Diameter <-
  #seperate out the initial sizes of urchins on day 1, when desired conditions were met.
  results %>% 
    select("Day", "Diam1", "Diam2", "TankID", "Temperature", "pH", "Treatment") %>%
    filter(TankID != "8t" ) %>% 
      #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
    filter(Day == "126") %>% 
      #sizes on day 1  
    select("TankID", "Diameter", "Temperature", "pH", "Treatment")


Day126Mean <- mean(Day126Diameter$Diameter)
Day126SE <- sd(Day126Diameter$Diameter)/sqrt(23)

Day126Mean
Day126SE
v. Table 3: Urchin diameter means
Day Mean Urchin Diameter (mm) se
-24 7.53 0.15
-13 10.38 0.23
1 16.03 0.36
126 70.52 1.43

b. Results of linear mixed model effect on diameter (day 1)

Day1Mod <- lmer(Diameter ~ Temperature*pH + (1|TankID), data=Day1Diameter)
anova(Day1Mod, type = 2)
## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Temperature    0.04177 0.04177     1    19  0.1922 0.66606  
## pH             0.38094 0.38094     1    19  1.7525 0.20128  
## Temperature:pH 0.89969 0.89969     1    19  4.1389 0.05611 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions of diameter model (Day 1)

#check assumptions
# 1. Normality of residuals
qqPlot(residuals(Day1Mod)) #normal

shapiro.test(residuals(Day1Mod)) #passes
hist(residuals(Day1Mod)) #normal

#check assumptions
# 2. Equal variances
leveneTest(residuals(Day1Mod)~InitialDiameter$Treatment) #Passes
plot(fitted(Day1Mod),resid(Day1Mod,type="pearson"),col="blue")

3. Response Analysis

a. Growth

i. Growth Analysis

Assemble data and calculate means using day 1 as initial diameter and calculate means

### GROWTH MODEL 1: using day 1 as the first official day of experiment when desired environmental conditions were met (24 days after getting urchins. In that 24 days, they were acclimated and then conditions ramped up)

Initial <-
  #seperate out the initial sizes of urchins on day 1, when desired conditions were met to standardize growth to this measurement.
  results %>% 
    select("Day", "Treatment", "Temperature", "pH", "Diam1", "Diam2", "TankID") %>%
    filter(TankID != "8t" ) %>% 
      #remove urchin number 8 which died halfway through
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column in order to make number of levels of each grouping factor < number of observations (as per error message previously recieved)
    filter(Day == "1") %>% 
      #sizes on day 1  
    select("Temperature", "pH", "Diameter", "TankID")

  
Final <- 
  #seperate out the final sizes of urchins on day 126
   results %>% 
    select("Day", "Treatment", "Temperature", "pH", "Diam1", "Diam2", "TankID") %>% 
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
      #gather Diam1 and Diam2 into single column
    filter(Day == "126") %>% 
      #filter for sizes on day 126
    select("Temperature", "pH", "Diameter","TankID")


Growth <- 
  #create column that is growth
  bind_cols(Initial, Final) %>% 
    #binds initial (Diameter) and final (Diameter1) sizes to calculate growth
  mutate(Growth = ((Diameter1 - Diameter)/Diameter*100)) %>% 
    #add column for growth calculation
  select("Growth")


resultsGrow <-
  #table of initial and final size data
  results %>% 
    select("Day", "Treatment", "Temperature", "pH", "TankID", "Diam1", "Diam2") %>% 
    drop_na(Diam1, Diam2) %>% 
    gather(key = "1or2", value = "Diameter", "Diam1", "Diam2") %>% 
    filter(Day == "126")
resultsGrow <- bind_cols(resultsGrow, Growth)
#combines table of initial and final with growth information

#Figure out means of each treatment group
resultsGrow %>% 
  dplyr::group_by(Temperature, pH) %>% 
  dplyr::summarise(mean = mean(Growth), s.e. = se(Growth)) %>% 

  kbl(caption = "Table 6: Growth means (%) +/- se after 126 days") %>% 
  kable_classic_2() %>% 
  kable_styling(c("striped", "hover"), full_width = F, html_font= "Times New Roman")
Table 6: Growth means (%) +/- se after 126 days
Temperature pH mean s.e.
ambient acidified 326.5950 16.290140
ambient ambient 323.1710 10.963618
heated acidified 352.0342 15.002253
heated ambient 376.3414 7.147267

Growth ANOVA results

#LMM for growth:
modGrow <-
  resultsGrow %>% 
  aov(Growth~ Temperature * pH, data=.)

summary(modGrow)
##                Df Sum Sq Mean Sq F value  Pr(>F)   
## Temperature     1  16750   16750   8.257 0.00634 **
## pH              1   1096    1096   0.540 0.46642   
## Temperature:pH  1   2197    2197   1.083 0.30396   
## Residuals      42  85202    2029                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions for growth model

## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modGrow)) 

## [1] 38 34
shapiro.test(residuals(modGrow)) #passes
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modGrow)
## W = 0.98242, p-value = 0.7064
hist(residuals(modGrow)) #looks normal

## ASSSUMPTIONS
# 2. Equal variances
leveneTest(residuals(modGrow)~resultsGrow$Treatment) #Passes
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  3  2.3134 0.08975 .
##       42                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fitted(modGrow),resid(modGrow,type="pearson"),col="blue")

Power analysis for test

pwr.anova.test(k=4, n= 5.75, f=0.7, sig.level = 0.05) #balanced one way anova
## 
##      Balanced one-way analysis of variance power calculation 
## 
##               k = 4
##               n = 5.75
##               f = 0.7
##       sig.level = 0.05
##           power = 0.7202639
## 
## NOTE: n is number in each group
#k=Number of groups - 4 treatments
#n=Number of observations (per group) - 5.75 to account for 1 dead
#f=Effect size - 0.5 based on some readings, but not sure if that's right
#sig.level=Significance level (Type I error probability)
#power=Power of test (1 minus Type II error probability)


pwr.f2.test(u=1, v=42, f=0.7, sig.level = 0.05) #general linear model
## 
##      Multiple regression power calculation 
## 
##               u = 1
##               v = 42
##              f2 = 0.7
##       sig.level = 0.05
##           power = 0.9997292
#u  = degrees of freedom for numerator
#v  = degrees of freedom for denominator
#f2 = effect size
#sig.level  = Significance level (Type I error probability)
#power  = Power of test (1 minus Type II error probability)

ii. Growth Interaction Plot

growthsummary<-resultsGrow %>% 
  #filter(Day==126) %>%
  select("pH", "Temperature","Growth") %>%
  drop_na(Growth) %>% 
  mutate(meanPct = Growth) %>%
  group_by(pH, Temperature) %>%
  summarise(mean = mean(meanPct), 
            N = length(meanPct),
            se = std.error(meanPct))
   
  
growthplot<-ggplot(data=growthsummary, aes(x=Temperature, y=mean, color = pH, group=interaction(Temperature,pH))) +
    geom_point(size = 3,  position=position_dodge(0.1)) +
    geom_line(aes(group=pH),  position=position_dodge(0.1)) +
    xlab("Temperature Treatment") +
    ylab("Total Growth (%)")+
    geom_errorbar(aes(ymin = mean-se, ymax = mean+se, width = 0.1), position=position_dodge(0.1)) +
    theme_bw() + 
    scale_x_discrete(limits=c("ambient", "heated"), 
                     label=c("Ambient", "High"))+
    scale_colour_manual(name = "pH Treatment",
                        values = c("skyblue3", "gray40"),
                        labels = c("Acidified", "Ambient"),
                        guide = guide_legend(reverse = TRUE)) +
    geom_text(x=1.5, y=125, label="p(pH)=0.466", size=6, color="darkgray") + 
    geom_text(x=1.5, y=100, label="p(Temperature)=0.006", size=6, color="black") + 
    geom_text(x=1.5, y=75, label="p(Interaction)=0.303", size=6, color="darkgray") + 
    ylim(0,450)+
    theme_classic()+
    theme(legend.position = "right", 
          legend.background = element_blank(),
          axis.title = element_text(size = 16, face="bold"),
          axis.text = element_text(size = 14), 
          plot.title = element_text(size=16, face="bold"),
          legend.text = element_text(size=14), 
          legend.title = element_text(size=16, face="bold"));growthplot

ggsave(filename="Figures/20200611/GrowthFig.png", plot=growthplot, dpi=300, width=8, height=6, units="in")

b. Growth over time analysis

resultsGrowTime <-  #create table of growth in treatments
  results %>% 
  select("Day", "TankID", "Treatment","Temperature","pH","Growth") %>%
  drop_na(Growth) %>% 
  group_by(Day, Treatment) %>%
  mutate(se = (std.error(Growth))*100) %>%
  mutate(meanPct = (mean(Growth)*100))

modelGrowTime <- lmer(Growth~Temperature*pH + (1|TankID), data=resultsGrowTime)
summary(modelGrowTime) #generate results
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: Growth ~ Temperature * pH + (1 | TankID)
##    Data: resultsGrowTime
## 
## REML criterion at convergence: 1181.4
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.66228 -0.87498  0.05818  0.82380  2.19934 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  TankID   (Intercept) 0.00     0.000   
##  Residual             1.27     1.127   
## Number of obs: 382, groups:  TankID, 24
## 
## Fixed effects:
##                              Estimate Std. Error        df t value
## (Intercept)                   1.77781    0.11504 378.00000  15.454
## Temperatureheated             0.09323    0.16269 378.00000   0.573
## pHambient                    -0.03635    0.16269 378.00000  -0.223
## Temperatureheated:pHambient   0.03893    0.23069 378.00000   0.169
##                             Pr(>|t|)    
## (Intercept)                   <2e-16 ***
## Temperatureheated              0.567    
## pHambient                      0.823    
## Temperatureheated:pHambient    0.866    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.707              
## pHambient   -0.707  0.500       
## Tmprtrhtd:H  0.499 -0.705 -0.705
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(modelGrowTime, type=2)
## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## Temperature    1.21056 1.21056     1   378  0.9529 0.3296
## pH             0.02757 0.02757     1   378  0.0217 0.8830
## Temperature:pH 0.03618 0.03618     1   378  0.0285 0.8661

i. Growth over time plot across all treatments

#### Growth Across Treatments
results %>% 
  select("Day", "Treatment","Temperature","pH","Growth") %>%
  drop_na(Growth) %>% 
  group_by(Day, Treatment) %>%
  mutate(se = (std.error(Growth))*100) %>%
  mutate(meanPct = (mean(Growth)*100)) %>% 

  
ggplot(aes(Day, meanPct, shape = Treatment, color = Treatment)) +
    geom_point(size = 2.5) +
    geom_line() +
    scale_x_continuous(name="Time (Day)", breaks = seq(0,130, by = 10)) +
    scale_y_continuous(name="Growth (%)", breaks = seq(0,400, by = 25)) +
    geom_errorbar(aes(ymin = meanPct-se, ymax = meanPct+se, width = 1.5)) +
    theme_classic() + 
    scale_shape_discrete(name = NULL,
                         labels = c("Amb T, Amb pH", "Amb T, Low pH", "High T, Amb pH", "High T, Low pH")) +
    scale_colour_manual(name = NULL,
                        values = c("gray40", "skyblue3", "firebrick3","mediumpurple3"),
                        labels = c("Amb T, Amb pH", "Amb T, Low pH", "High T, Amb pH", "High T, Low pH")) +
    theme(legend.margin = margin(0), 
          legend.position = c(.1, .98), 
          legend.justification = c("left", "top"),
          legend.background = element_blank())

ggsave(filename="Figures/20200611/GrowthTreatments.png", dpi=300, width=8, height=6, units="in")

ii. Growth over time plot by temperature

results %>% 
  select("Day","Temperature","Growth") %>%
  drop_na(Growth) %>% 
  group_by(Day, Temperature) %>%
  mutate(se = (std.error(Growth))*100) %>%
  mutate(meanPct = (mean(Growth)*100)) %>% 

  
  ggplot(aes(Day, meanPct, shape = Temperature, color = Temperature)) +
    geom_point(size = 2.5) +
    geom_line() +
    scale_x_continuous(name="Time (Day)", breaks = seq(0,130, by = 10)) +
    scale_y_continuous(name="Growth (%)", breaks = seq(0,400, by = 25)) +
    geom_errorbar(aes(ymin = meanPct-se, ymax = meanPct+se, width = 1.5)) +
    theme_classic() +
    scale_colour_manual(name = NULL, 
                        values = c("gray40", "firebrick3"), 
                        labels = c("Ambient Temperture", "High Temperature")) +
    scale_shape_discrete(name = NULL, 
                       labels = c("Ambient Temperture", "High Temperature")) +
    theme(legend.margin = margin(0), 
          legend.position = c(.03, .93), 
          legend.justification = c("left", "top")) + 
    annotate("text", x = 0, y = 375, label = "(a)") +
    annotate("text", x = 130, y = 359, label = "*", size = 7)

iii. Growth over time plot by pH

results %>% 
  select("Day","pH","Growth") %>%
  drop_na(Growth) %>% 
  group_by(Day, pH) %>%
  mutate(se = (std.error(Growth))*100) %>%
  mutate(meanPct = (mean(Growth)*100)) %>% 

  ggplot(aes(Day, meanPct, shape = pH, color = pH)) +
    geom_point(size = 2.5) +
    geom_line() +
    scale_x_continuous(name="Time (Day)", breaks = seq(0,130, by = 10)) +
    scale_y_continuous(name="Growth (%)", breaks = seq(0,400, by = 25)) +
    geom_errorbar(aes(ymin = meanPct-se, ymax = meanPct+se, width = 1.5)) +
    theme_classic() +
    scale_colour_manual(name = NULL, 
                        values = c("skyblue3", "gray40"), 
                        labels = c("Low pH", "Ambient pH")) +
    scale_shape_discrete(name = NULL,
                         labels = c("Low pH", "Ambient pH")) +
    theme(legend.margin = margin(0), 
          legend.position = c(.03, .93), 
          legend.justification = c("left", "top")) + 
    annotate("text", x = 0, y = 370, label = "(b)")

b. Calcification Ratio

i. Calcification Ratio Analysis at the tips

Assemble data for chosen spine tips (distal end)

### Model 1: calcification at the tips of spines

resultsTip <- 
  #create data using only the SEM images of the tips of spines.
  results %>% 
  select ("Treatment", "Temperature", "pH", "PartOfSpine", "Chosen", "RatioSEM", "TankID") %>% 
    #selects for desired variables. PartOfSpine is tip or base, Chosen indicates which side of the image (arbitrary) is selected for (aka has no dust)
  drop_na(RatioSEM) %>% 
  filter(Chosen == "yes", PartOfSpine == "Tip")
    #filters only for tips and chosen side. No dusty images to be analyzed.

Results for spine tips with linear mixed model

#LMM for calcification at the tips of spines 
modTip <-
  resultsTip %>% 
  lmer(RatioSEM ~ Temperature * pH + (1|TankID), data = .)

summary(modTip)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RatioSEM ~ Temperature * pH + (1 | TankID)
##    Data: .
## 
## REML criterion at convergence: 64.1
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.76413 -0.74507 -0.02052  0.63045  2.07630 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  TankID   (Intercept) 0.0000   0.0000  
##  Residual             0.1356   0.3682  
## Number of obs: 67, groups:  TankID, 23
## 
## Fixed effects:
##                             Estimate Std. Error       df t value Pr(>|t|)
## (Intercept)                  1.69118    0.08930 63.00000  18.938   <2e-16
## Temperatureheated           -0.06662    0.12453 63.00000  -0.535   0.5945
## pHambient                   -0.07462    0.12453 63.00000  -0.599   0.5512
## Temperatureheated:pHambient  0.30557    0.18089 63.00000   1.689   0.0961
##                                
## (Intercept)                 ***
## Temperatureheated              
## pHambient                      
## Temperatureheated:pHambient .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.717              
## pHambient   -0.717  0.514       
## Tmprtrhtd:H  0.494 -0.688 -0.688
## convergence code: 0
## boundary (singular) fit: see ?isSingular
anova(modTip, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type II Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF DenDF F value  Pr(>F)  
## Temperature    0.10158 0.10158     1    63  0.7492 0.39000  
## pH             0.08185 0.08185     1    63  0.6038 0.44006  
## Temperature:pH 0.38685 0.38685     1    63  2.8534 0.09612 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions for calcification at tip model.

## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modTip)) #Normal

## [1] 38  8
shapiro.test(residuals(modTip)) #Pass
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modTip)
## W = 0.97737, p-value = 0.2605
hist(residuals(modTip))

# 2. Equal variances
leveneTest(residuals(modTip)~resultsTip$Treatment) #Pass
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3  0.4545  0.715
##       63
plot(fitted(modTip),resid(modTip,type="pearson"),col="blue") 

ii. Calcification Ratio Analysis at the base

Assemble data for chosen spine bases (proximal end)

### Model 2: calcification in the base of the spines

resultsBase <-
   #create data using only the SEM images of the base of spines.
  results %>% 
  select ("Treatment", "Temperature", "pH", "PartOfSpine", "Chosen", "RatioSEM", "TankID") %>% 
    #selects for desired variables. PartOfSpine is tip or base, Chosen indicates which side of the image (arbitrary) is selected for (aka has no dust)
  drop_na(RatioSEM) %>% 
  filter(Chosen == "yes", PartOfSpine == "Base")
    #filters only for tips and chosen side. No dusty images to be analyzed.

Results for spine bases with linear mixed model

#LMM for calcification at the tips of spines 
modBase <- 
  resultsBase %>% 
  lmer(RatioSEM~ Temperature * pH + (1|TankID), data=.)

summary(modBase) 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: RatioSEM ~ Temperature * pH + (1 | TankID)
##    Data: .
## 
## REML criterion at convergence: 98.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.5535 -0.6315 -0.1529  0.6918  1.9937 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  TankID   (Intercept) 0.1303   0.3610  
##  Residual             0.1597   0.3996  
## Number of obs: 68, groups:  TankID, 23
## 
## Fixed effects:
##                             Estimate Std. Error      df t value Pr(>|t|)
## (Intercept)                   2.1440     0.1749 18.8077  12.257 2.06e-10
## Temperatureheated             0.3566     0.2487 19.1707   1.434   0.1677
## pHambient                     0.6524     0.2474 18.8077   2.637   0.0163
## Temperatureheated:pHambient  -0.2224     0.3594 18.9804  -0.619   0.5434
##                                
## (Intercept)                 ***
## Temperatureheated              
## pHambient                   *  
## Temperatureheated:pHambient    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Tmprtr pHmbnt
## Tempertrhtd -0.703              
## pHambient   -0.707  0.497       
## Tmprtrhtd:H  0.487 -0.692 -0.688
anova(modBase, type="II") #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.
## Type II Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Temperature    0.30743 0.30743     1 18.992  1.9251 0.181362   
## pH             1.48873 1.48873     1 18.960  9.3223 0.006553 **
## Temperature:pH 0.06114 0.06114     1 18.980  0.3829 0.543424   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions for calcification at base model.

## ASSSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modBase)) #Normal

## [1] 33 37
shapiro.test(residuals(modBase)) # Pass
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modBase)
## W = 0.97304, p-value = 0.1472
hist(residuals(modBase))

# 2. Equal variances
leveneTest(residuals(modBase)~resultsBase$Treatment) #Pass
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value Pr(>F)
## group  3   1.068  0.369
##       64
plot(fitted(modBase),resid(modBase,type="pearson"),col="blue")

iii. Figure 4: Calcification Ratio Interaction Plot at tips and bases

Make Interaction plots for calcification at the base and the tip

#### Tip and base calcification ratios

#base
basesummary<-results %>% 
  select("Temperature","pH","RatioSEM", "PartOfSpine") %>%
  filter(PartOfSpine=="Base")%>%
  drop_na(RatioSEM) %>% 
  group_by(pH, Temperature, PartOfSpine) %>%
  summarize(mean = mean(RatioSEM), se = std.error(RatioSEM), N = length(RatioSEM))
  
baseplot<- ggplot(data=basesummary, aes(x=Temperature, y=mean, color = pH, group=interaction(Temperature,pH))) +
    geom_point(size = 3,  position=position_dodge(0.1)) +
    geom_line(aes(group=pH),  position=position_dodge(0.1)) +
    xlab("Temperature Treatment") +
    ylab("Proximal Calcification Ratio")+
    ylim(1,3.5)+
    geom_errorbar(aes(ymin = mean-se, ymax = mean+se, width = 0.1), position=position_dodge(0.1)) +
    theme_bw() + 
    scale_x_discrete(limits=c("ambient", "heated"), 
                     label=c("Ambient", "High"))+
    scale_colour_manual(name = "pH Treatment",
                        values = c("skyblue3", "gray40"),
                        labels = c("Acidified", "Ambient"),
                        guide = guide_legend(reverse = TRUE)) +
    geom_text(x=1.5, y=1.6, label="p(pH)=0.002", size=6, color="black") + 
    geom_text(x=1.5, y=1.4, label="p(Temperature)=0.16", size=6, color="darkgray") + 
    geom_text(x=1.5, y=1.2, label="p(Interaction)=0.54", size=6, color="darkgray") + 
    theme_classic()+
    theme(legend.position = "right", 
          legend.background = element_blank(), 
          axis.title = element_text(size = 16, face="bold"),
          axis.text = element_text(size = 14), 
          plot.title = element_text(size=16, face="bold"),
          legend.text = element_text(size=14), 
          legend.title = element_text(size=16, face="bold"));baseplot

#ggsave("CalcRatiobase.png", plot=baseplot, path = "/Users/emilysesno/Desktop/R_Analysis/R_Analysis/output",  width = 7, height = 5)  

ggsave(filename="Figures/20200611/BaseCalcificationFig.png", plot=baseplot, dpi=300, width=8, height=6, units="in")

#tip
tipsummary<-results %>% 
  select("Temperature","pH","RatioSEM", "PartOfSpine") %>%
  filter(PartOfSpine=="Tip")%>%
  drop_na(RatioSEM) %>% 
  group_by(pH, Temperature, PartOfSpine) %>%
  summarize(mean = mean(RatioSEM), se = std.error(RatioSEM), N = length(RatioSEM))

tipplot<-ggplot(data=tipsummary, aes(x=Temperature, y=mean, color = pH, group=interaction(Temperature,pH))) +
    geom_point(size = 3,  position=position_dodge(0.1)) +
    geom_line(aes(group=pH),  position=position_dodge(0.1)) +
    xlab("Temperature Treatment") +
    ylab("Distal Calcification Ratio")+
    ylim(1,3.5)+
    geom_errorbar(aes(ymin = mean-se, ymax = mean+se, width = 0.1), position=position_dodge(0.1)) +
    theme_bw() + 
    scale_x_discrete(limits=c("ambient", "heated"), 
                     label=c("Ambient", "High"))+
    scale_colour_manual(name = "pH Treatment",
                        values = c("skyblue3", "gray40"),
                        labels = c("Acidified", "Ambient"),
                        guide = guide_legend(reverse = TRUE)) +
    geom_text(x=1.5, y=3.0, label="p(pH)=0.44", size=6, color="darkgray") + 
    geom_text(x=1.5, y=2.8, label="p(Temperature)=0.39", size=6, color="darkgray") + 
    geom_text(x=1.5, y=2.6, label="p(Interaction)=0.09", size=6, color="darkgray") + 
    theme_classic()+
    theme(legend.position = "none", #removed the legend for the tip, so when we align them horizontally we only have one legend. But need to adjust the size so it is the same as base
          legend.background = element_blank(), 
          axis.title = element_text(size = 16, face="bold"),
          axis.text = element_text(size = 14), 
          plot.title = element_text(size=16, face="bold"),
          legend.text = element_text(size=14), 
          legend.title = element_text(size=16, face="bold"));tipplot

#ggsave(filename="Figures/20200611/TipCalcificationFig.png", plot=tipplot, dpi=300, width=8, height=6, units="in")

Combine both interaction plots at base and tip to be in horizontal line

c. Spine Loss

i. Spine Loss Analysis:

Assemble data for dropped spines

## Number of loose spines counted at the bottom of tanks at the end of the experiment. These were either broken or shed. Note: these spines were not collected, just counted, so can not determine whether they were broken or shed.

resultsDropped <-
  #create data using only the needed variables.
  results %>% 
  select ("Treatment", "Temperature", "pH", "SpineCount", "TankID") %>% 
  drop_na(SpineCount) %>% 
  group_by(pH) %>% 
  mutate(mean = mean(SpineCount),
         se = std.error(SpineCount))

Analyse dropped spines with linear mixed effect model.

##LM of dropped spines - can't make it LMM because error: number of levels of each grouping factor must be < number of observations - only have one count for each TankID.
modDropped <- 
  resultsDropped %>% 
  lm(SpineCount~ Temperature * pH, data=.)

summary(modDropped)
## 
## Call:
## lm(formula = SpineCount ~ Temperature * pH, data = .)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -12.0000  -6.0833  -0.1667   0.4000  20.0000 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                   21.167      3.640   5.814 1.33e-05 ***
## Temperatureheated             -1.167      5.148  -0.227  0.82315    
## pHambient                    -21.000      5.148  -4.079  0.00064 ***
## Temperatureheated:pHambient    1.600      7.461   0.214  0.83248    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.917 on 19 degrees of freedom
## Multiple R-squared:  0.6088, Adjusted R-squared:  0.547 
## F-statistic: 9.855 on 3 and 19 DF,  p-value: 0.0003899
anova(modDropped) #anova{stats}.Compute analysis of variance (or deviance) tables for one or more fitted model objects. Produces Type III Analysis of Variance Table with Satterthwaite's method.
## Analysis of Variance Table
## 
## Response: SpineCount
##                Df  Sum Sq Mean Sq F value    Pr(>F)    
## Temperature     1    1.52    1.52  0.0192    0.8914    
## pH              1 2345.78 2345.78 29.4995 3.062e-05 ***
## Temperature:pH  1    3.66    3.66  0.0460    0.8325    
## Residuals      19 1510.87   79.52                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Check assumptions for model for dropped spines

##ASSUMPTIONS
# 1. Normality of residuals
qqPlot(residuals(modDropped)) #not normal

## [1] 2 4
shapiro.test(residuals(modDropped)) #faaaiiilll
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modDropped)
## W = 0.86706, p-value = 0.005636
hist(residuals(modDropped))

# 2. Equal variances
leveneTest(residuals(modDropped)~resultsDropped$Treatment) # barely pass
## Levene's Test for Homogeneity of Variance (center = median)
##       Df F value  Pr(>F)  
## group  3  2.8487 0.06483 .
##       19                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(fitted(modDropped),resid(modDropped,type="pearson"),col="blue")

Analysis for dropped spines does not meet assumptions. Attempt a square root transformation.

##DOES NOT MEET ASSUMPTIONS so ATTEMPT AT Transformation:
resultsDropped$tdata<-(resultsDropped$SpineCount)^1/2
  # transformation
modDropped1 <- lm(tdata~ Temperature * pH, data=resultsDropped)
  # run lm model of transformed data

###ASSUMPTIONS for Transformation
# 1. Normality of residuals
qqPlot(residuals(modDropped1)) #not normal

## [1] 2 4
shapiro.test(residuals(modDropped1)) #fail
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(modDropped1)
## W = 0.86706, p-value = 0.005636
hist(residuals(modDropped1))

# 2. equal variances
bartlett.test(residuals(modDropped1)~resultsDropped$Treatment) #fail
## 
##  Bartlett test of homogeneity of variances
## 
## data:  residuals(modDropped1) by resultsDropped$Treatment
## Bartlett's K-squared = 43.282, df = 3, p-value = 2.144e-09
plot(fitted(modDropped1),resid(modDropped1,type="pearson"),col="blue")

summary(modDropped1)
## 
## Call:
## lm(formula = tdata ~ Temperature * pH, data = resultsDropped)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.0000 -3.0417 -0.0833  0.2000 10.0000 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  10.5833     1.8202   5.814 1.33e-05 ***
## Temperatureheated            -0.5833     2.5742  -0.227  0.82315    
## pHambient                   -10.5000     2.5742  -4.079  0.00064 ***
## Temperatureheated:pHambient   0.8000     3.7304   0.214  0.83248    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.459 on 19 degrees of freedom
## Multiple R-squared:  0.6088, Adjusted R-squared:  0.547 
## F-statistic: 9.855 on 3 and 19 DF,  p-value: 0.0003899
Anova(modDropped1, type=2) #test for significance of model
## Anova Table (Type II tests)
## 
## Response: tdata
##                Sum Sq Df F value    Pr(>F)    
## Temperature      0.23  1  0.0118    0.9146    
## pH             586.44  1 29.4995 3.062e-05 ***
## Temperature:pH   0.91  1  0.0460    0.8325    
## Residuals      377.72 19                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Transfomraiton (of any usual kind) of dropped spines was not successful, so we use a non parametric Kruskal Wallis Test.

### Transformations not work, so use non parametric:
kruskal.test(resultsDropped$SpineCount~resultsDropped$Treatment)# SIG*** (p=0.0005563)

kruskal.test(resultsDropped$SpineCount~resultsDropped$Temperature)# 
kruskal.test(resultsDropped$SpineCount~resultsDropped$pH)#

Conduct Dunn Post Hoc Test for dropped spines

## POST HOC Pairwise
dunn <- 
  resultsDropped %>% 
  dunnTest(SpineCount ~ Treatment,
           method = "holm", kw = TRUE, data = .)
dunnph <- dunn$res
cldList(P.adj ~ Comparison, data = dunnph, threshold = 0.05)
##   Group Letter MonoLetter
## 1   AMB      a        a  
## 2    OA      b         b 
## 3     T     ac        a c
## 4  T/OA     bc         bc

ii. Figure 5: Spine Loss Interaction Plot

droppedsummary<-results %>% 
  select("Temperature","pH","SpineCount") %>%
  drop_na(SpineCount) %>% 
  group_by(pH, Temperature) %>%
  summarize(mean = mean(SpineCount), se = std.error(SpineCount), N = length(SpineCount))

  dropplot<-ggplot(data=droppedsummary, aes(x=Temperature, y=mean, color = pH, group=interaction(Temperature,pH))) +
    geom_point(size = 3,  position=position_dodge(0.1)) +
    geom_line(aes(group=pH),  position=position_dodge(0.1)) +
    xlab("Temperature Treatment") +
    ylab("Dropped Spines")+
    geom_errorbar(aes(ymin = mean-se, ymax = mean+se, width = 0.1), position=position_dodge(0.1)) +
    ylim(0,30) +
    theme_bw() + 
    scale_x_discrete(limits=c("ambient", "heated"), 
                     label=c("Ambient", "High"))+
    scale_colour_manual(name = "pH Treatment",
                        values = c("skyblue3", "gray40"),
                        labels = c("Acidified", "Ambient"),
                        guide = guide_legend(reverse = TRUE)) +
    #geom_text(x=1.5, y=11, label="p(pH)<0.001", size=6, color="black") + 
    #geom_text(x=1.5, y=9, label="p(Temperature)=0.91", size=6, color="darkgray") + 
    #geom_text(x=1.5, y=7, label="p(Interaction)=0.83", size=6, color="darkgray") +
    theme_classic()+
    theme(legend.position = "right", 
          legend.background = element_blank(), 
          axis.title = element_text(size = 16, face="bold"),
          axis.text = element_text(size = 14), 
          plot.title = element_text(size=16, face="bold"),
          legend.text = element_text(size=14), 
          legend.title = element_text(size=16, face="bold"));dropplot

#ggsave(filename="Figures/20200611/DroppedSpinesFig.png", plot=dropplot, dpi=300, width=8, height=6, units="in")